#### Created by Yu Luo####
#Created March, 20, 2021
#Last updated April, 24, 2023
#Purpose: meta-analysis of nudge and sludge interventions by cognitive processes
#Code source:
  #Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019). 
  #Doing Meta-Analysis in R: A Hands-on Guide. 
  #https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/.

library(dplyr)
library(esc)
# please change the cells from Kettles et al. from general to numbers in the csv to avoid errors
d <- read.csv("rawdata_included_only.csv", encoding='utf-8')
d <- filter(d, Included==1)


metaD <- data.frame()

#Convert the effect sizes to Cohen's d
for (i in 1:nrow(d)){
  
  if (d$Function[i]=="esc_mean_sd"){
    currentStudy <- esc_mean_sd(grp1m = d$grp1m[i], grp1sd = d$grp1sd[i], grp1n = d$grp1n[i], 
                                grp2m = d$grp2m[i], grp2sd = d$grp2sd[i], grp2n = d$grp2n[i], 
                                es.type = "d", study = d$Author[i])
    
  }else if (d$Function[i]=="esc_2x2"){
    currentStudy <- esc_2x2(grp1yes = d$grp1yes[i], grp1no = d$grp1no[i], 
                            grp2yes = d$grp2yes[i], grp2no = d$grp2no[i], 
                            es.type = "d", study = d$Author[i])

  }else if (d$Function[i]=="esc_bin_prop"){
    currentStudy <- esc_bin_prop(prop1event = d$prop1event[i], grp1n = d$grp1n.1[i], 
                                 prop2event = d$prop2event[i], grp2n = d$grp2n.1[i], 
                                 es.type = "d", study = d$Author[i])

  }else if (d$Function[i]=="esc_B"){
    currentStudy <- esc_B( d$b[i], sdy = d$sdy[i], 
                          grp1n =  d$grp1n.3[i],
                          grp2n =  d$grp2n.3[i], 
                          es.type = "d", study = d$Author[i])
    
  }else if (d$Function[i]=="esc_beta"){
    currentStudy <- esc_beta( d$B[i], sdy = d$sdy.1[i], 
                           grp1n =  d$grp1n.6[i],
                           grp2n =  d$grp2n.6[i], 
                           es.type = "d", study = d$Author[i])
  
  }else if (d$Function[i]=="esc_mean_se"){
    currentStudy <- esc_mean_se(grp1m = d$grp1m.1[i], grp1se = d$grp1se[i], grp1n = d$grp1n.4[i], 
                                grp2m = d$grp2m.1[i], grp2se = d$grp2se[i], grp2n = d$grp2n.4[i], 
                                es.type = "d", study = d$Author[i])
    
  }else if (d$Function[i]=="esc_f"){
    currentStudy <- esc_f(f=d$F[i], grp1n=d$grp1n.2[i], grp2n=d$grp2n.2[i], 
                          es.type = "d", study = d$Author[i])

  }else if (d$Function[i]=="esc_t"){
    currentStudy <- esc_t(t = d$t[i], p=d$p[i], grp1n = d$grp1n.5[i], grp2n = d$grp2n.5[i], 
                          es.type = "d", study = d$Author[i])
  
  }else if (d$Function[i]=="convert_or2d"){
    currentStudy <- convert_or2d(or = d$or[i], 
                                 se = d$se.1[i], 
                                 totaln=d$totaln.3[i], 
                                 study = d$Author[i])
  
  }else if (d$Function[i]=="esc_mean_gain"){
    currentStudy <-esc_mean_gain(pre1mean = d$pre1mean[i], pre1sd = d$pre1sd[i], post1mean = d$post1mean[i],
                                 post1sd = d$post1sd[i], grp1n = d$grp1n.7[i], pre2mean = d$pre2mean[i], pre2sd = d$pre2sd[i],
                                 post2mean = d$post2mean[i], post2sd = d$post2sd[i], grp2n = d$grp2n.7[i], study = d$Author[i])
  }
  
  currentStudy <- data.frame(currentStudy)

  #Reverse coding for studies that aimed to inhibit a behavior (e.g., reduce water consumption) 
  if(d$Reversed[i]=="y"){
    currentStudy$reversed <- 1
    currentStudy$convertedES <- currentStudy$es*(-1)
  }else{
    currentStudy$reversed <- 0
    currentStudy$convertedES <- currentStudy$es
  }
  

  currentStudy$Cognitive_process <- d$Cognitive.process[i]
  currentStudy$Intervention <- d$Intervention[i]
  currentStudy$Domain <- d$Domain[i]
  currentStudy$DurationConverted <- d$DurationConverted[i]
  currentStudy$Perspective <- d$Perspective[i]
  currentStudy$Nudge_type <- d$Type[i]
  currentStudy$Behavior <- d$Behavioral.measure[i]
  currentStudy$Duration <- d$Duration[i]
  currentStudy$Combined <- d$Combined[i]
  currentStudy$Conclusion <- d$Findings[i]
  currentStudy$Year <- d$Year[i]
  
  
  
  metaD <- rbind(metaD, currentStudy)
  
}

# write.csv(metaD, "meta-analysisData.csv", row.names = F)
write.csv(metaD, "meta-analysis_output.csv", row.names = F)


### descriptive statistics)
#number of studies per year
table(metaD$Nudge_type)
num_by_year <- as.data.frame(table(metaD$Year))

library(ggplot2)
#plot the number of studies per year
year_graph <- ggplot(num_by_year, aes(x = Var1, y = Freq)) + 
  geom_bar(stat = "identity",  fill = "#BFBFBF", width = 0.6) +
  ylab("Number of articles") +
  xlab("")+
  theme_classic()+
  theme(text=element_text(size=20,  family="serif"),
        axis.text.x = element_text(color = "black"),
        axis.text.y = element_text(color = "black"))+
  scale_y_continuous(expand = c(0, 0), limits = c(0, 40),breaks = seq(0, 40, 5))

ggsave(file=paste("figure/year_plot", ".png", sep = "", collapse = NULL), plot=year_graph, width=8, height=4)


# sample size
sum(metaD$sample.size)


####Analyze the effect sizes (d)####
library(meta)
library(metafor)

newMeta <- read.csv("meta-analysis_output.csv")


#get the mean and se
library(dplyr)
dataoutput <- newMeta %>%
  group_by(Domain, Cognitive_process) %>%
  summarise(
    mean = mean(convertedES, na.rm = T),
    sd = sd(convertedES, na.rm = T),
    n = n(),
    se = sd / sqrt(n)
  )

write.csv(dataoutput, "table/meanbyDomain&Cogproc.csv", row.names = F)


metaData <- metagen(TE=newMeta$convertedES,
                  seTE=newMeta$se,
                  data = newMeta,
                  studlab = paste(study),
                  comb.fixed = TRUE,
                  comb.random = TRUE,
                  method.tau = "SJ",
                  hakn = T,
                  prediction = TRUE,
                  sm = "SMD")
metaData
summary(metaData)

#function to build the output file
getTable<- function(variableNeeded, outputFileName, ordering){
  #by cognitive process and nudge/sludge
  cogproc.subgroup<-update.meta(metaData, 
                                 byvar=eval(parse(text = variableNeeded)), 
                                 comb.random = TRUE, 
                                 comb.fixed = FALSE)
  
  print(summary(cogproc.subgroup))
  
  # cogproc.subgroup
  # print(cogproc.subgroup, digits = 2)
  
  png(file=paste("plot/", outputFileName, ".png"),
      width=1250, height=4000)
  forest(cogproc.subgroup)
  dev.off()
  
  
  #build the table
  CI <- NA
  
  cogProcess <- data.frame(cogproc.subgroup$bylevs, 
                           cogproc.subgroup$k.w, 
                           CI,
                           cogproc.subgroup$TE.random.w, 
                           cogproc.subgroup$lower.random.w,
                           cogproc.subgroup$upper.random.w)
  
  cogProcess$CI <- paste(round(cogproc.subgroup$TE.random.w, 2), 
                         " [", round(cogproc.subgroup$lower.random.w, 2), ", ", 
                         round(cogproc.subgroup$upper.random.w, 2), "]", sep="")
  
  cogProcess <- cogProcess[order(match(cogProcess[[1]], ordering)), ]
  
  write.csv(cogProcess, paste("table/", outputFileName, ".csv", sep=""), row.names = FALSE)
  
  return(cogProcess)
}

#by nudge type
ordering <- c("Nudge", "Sludge")
getTable("Nudge_type", "byType", ordering)

#by cognitive process
ordering <- c("Attention", "Perception", "Memory", "Effort", "Intrinsic", "Extrinsic")
getTable("Cognitive_process", "byCogProcess", ordering)


#by cognitive process and type(nudge/sludge)
ordering <- c("Attention_Nudge", 
              "Attention_Sludge",
              "Perception_Nudge",
              "Perception_Sludge",
              "Memory_Nudge",
              "Memory_Sludge",
              "Effort_Nudge",
              "Effort_Sludge",
              "Intrinsic_Nudge",
              "Intrinsic_Sludge",
              "Extrinsic_Nudge",
              "Extrinsic_Sludge")
getTable("Combined", "byCogProcess&Type", ordering)


#by Domain
outputbyDom <- getTable("Domain", "byDomain", ordering)
outputbyDom <- outputbyDom[order(-outputbyDom$cogproc.subgroup.TE.random.w),]
write.csv(outputbyDom, "table/byDomain.csv", row.names = F)

#by intervention
outputbyInt <- getTable("Intervention", "intervention", ordering)
outputbyInt <- filter(outputbyInt, cogproc.subgroup.k.w >1)
outputbyInt <- outputbyInt[order(-outputbyInt$cogproc.subgroup.TE.random.w),]
write.csv(outputbyInt, "table/byIntervention.csv", row.names = F)

##### Meta regression #####

library(dplyr)
library(metafor)
library(multcomp)
regressionData <- read.csv("meta-analysis_output.csv")
regressionData$Cognitive_process  <- factor(regressionData$Cognitive_process , ordered = F)
regressionData$Cognitive_process  <- relevel(regressionData$Cognitive_process , ref = "Intrinsic")
ordering <- c("Intrinsic", "Effort", "Attention", "Extrinsic", "Perception", "Memory")
regressionData$Cognitive_process <- factor(regressionData$Cognitive_process, levels = ordering)

metaRegCog <- function(moderator, dataReg){
  model2 <- rma(yi = convertedES, 
                sei = se, 
                data = dataReg, 
                method = "ML", 
                mods = ~ eval(rlang::parse_expr(moderator)),
                test="knha")
  print(model2)
  
  print(summary(glht(model2, linfct=cbind(contrMat(c("Intrinsic"=1,
                                                     "Effort"=1,
                                                     "Attention"=1,
                                                     "Extrinsic"=1,
                                                     "Perception"=1,
                                                     "Memory"=1), type="Tukey"))), test=adjusted("none")))
  
  print(levels(regressionData$Cognitive_process))
}

#cognitive process
metaRegCog('Cognitive_process', regressionData)

#multivariate meta regression
metaRegSingle <- function(mod1, dataReg){
  model2 <- rma(yi = convertedES, 
                sei = se, 
                data = dataReg, 
                method = "ML", 
                mods = ~ eval(rlang::parse_expr(mod1)),
                test="knha")
  return(model2)
}

#nudge vs sludge
metaRegSingle('Nudge_type', regressionData)

#by domain:  cognitive process

#"environment"
environmentReg <- filter(regressionData, Domain=="environment")
metaRegCog('Cognitive_process', environmentReg)

#"finance"
financeReg <- filter(regressionData, Domain=="finance")
metaRegCog('Cognitive_process', financeReg)

#"education"
educationReg <- filter(regressionData, Domain=="education")
metaRegCog('Cognitive_process', educationReg) # multiple cognitive processes have no interventions

#"health" 
healthReg <- filter(regressionData, Domain=="health")
metaRegCog('Cognitive_process', healthReg)

#"policy"     
policyReg <- filter(regressionData, Domain=="policy")
metaRegCog('Cognitive_process', policyReg)

#by domain: nudge vs sludge
#"environment"
metaRegSingle('Nudge_type', environmentReg)

#"finance"
metaRegSingle('Nudge_type', financeReg)

#"education" - no sludge intervention

#"health" 
metaRegSingle('Nudge_type', healthReg)

#"policy"     
metaRegSingle('Nudge_type', policyReg)


#by Intervention
regressionDataMoreThan2 <-filter(regressionData, !Intervention %in% c("Alert", "Color", "Fine", "Font size", "Identifiable victim")) 
regressionDataMoreThan2$Intervention  <- factor(regressionDataMoreThan2$Intervention , ordered = F)
regressionDataMoreThan2$Intervention  <- relevel(regressionDataMoreThan2$Intervention , ref = "Social norm")
modelInt <- metaRegSingle('Intervention', regressionDataMoreThan2)
modelInt
levels(as.factor(regressionDataMoreThan2$Intervention))
output <- summary(glht(modelInt, linfct=cbind(contrMat(c(
   "Social norm"=1,              "Accessibility "=1,           "Active choice"=1,                               
   "Anchoring"=1,                "Appearance"=1,               "Availability"=1,                                
   "Commitment making"=1,        "Conditional incentives"=1,   "Convenience "=1,             "Default"=1,                 
   "Financial incentives"=1,                                         "Gain framing"=1,            
   "Goal setting"=1,             "Graphic"=1,                  "Highlighting"=1,                 
   "Implementation intention"=1, "Inconvenience"=1,            "Informational feedback"=1,   "Informational messaging"=1, 
   "Loss framing"=1,             "Motivational intervention"=1,"Non-financial incentives"=1, "Assortment size"=1,            
   "Priming"=1,                  "Reminder"=1,                 "Simplification"=1,           "Visibility"=1
), type="Tukey"))), test=adjusted("none"))

library(broom)
outputClean <- tidy(output)
outputClean <- subset(outputClean, select = -c(null.value, statistic))
outputClean$contrast <- ifelse(outputClean$estimate <0, 
                             gsub("(.*)-(.*)", "\\2-\\1", outputClean$contrast),
                             outputClean$contrast)
outputClean$estimate <- ifelse(outputClean$estimate <0, outputClean$estimate*-1, outputClean$estimate)

outputClean %>% 
            mutate(contrast = trimws(contrast))%>% 
            filter(p.value<.05) %>%
            mutate_if(is.numeric, round, 3) %>%
            arrange(desc(estimate)) %>%
            mutate(p.value = ifelse(p.value == 0, "<.001", p.value))%>%
            write.csv("table/table4TukeyComparision.csv", row.names = F)


#multivariate meta regression
metaRegCogType <- function(mod1, mod2, dataReg){
  model2 <- rma(yi = convertedES, 
                sei = se, 
                data = dataReg, 
                method = "ML", 
                mods = ~ eval(rlang::parse_expr(mod1))*eval(rlang::parse_expr(mod2)),
                test="knha")
  print(model2)
}

metaRegCogType('Cognitive_process', 'Nudge_type', regressionData)

#### Plot figures ####
library(ggpubr)
library(svglite)

ordering <- c("Effort", "Extrinsic", "Perception", "Attention", "Memory", "Intrinsic")
regressionData$Cognitive_process <- factor(regressionData$Cognitive_process, levels = ordering)
regressionData$Nudge_type <- factor(regressionData$Nudge_type, levels = c("Nudge", "Sludge"))

barPlotJitter <- function(iv, d, graphname){
  
  cog_count <- as.data.frame(table(d[[iv]]))
  cog_count$level_name <- paste0(cog_count$Var1, "\n(k=", cog_count$Freq, ")")
  levels(d[[iv]]) <- c(cog_count$level_name)

  graph <- ggbarplot(d, x = iv, y = c("convertedES"), 
                    merge = TRUE, add = c("mean_se"), 
                    ylab="Effect size (d)", 
                    xlab = "",
                    ylim = c(-2,3))+
          scale_x_discrete(drop = FALSE) +
          geom_jitter(color = "darkgray", width = 0.2, shape=16)+
          theme(text=element_text(size=20,  family="serif"))
  
  
  print(graph)
  
  if (iv=="Nudge_type"){
    ggsave(file=paste("figure/", graphname, ".svg", sep = "", collapse = NULL), plot=graph, width=3, height=5)
    
  }else{
    ggsave(file=paste("figure/", graphname, ".svg", sep = "", collapse = NULL), plot=graph, width=9, height=5)
    
  }
  
}

barPlotJitter("Cognitive_process", regressionData, "overall")

barPlotJitter('Nudge_type', regressionData, "nudge")

library(dplyr)
#"education"
educationReg <- filter(regressionData, Domain=="education")
educationReg$Nudge_type <- factor(educationReg$Nudge_type, levels = c("Nudge", "Sludge"))

barPlotJitter('Cognitive_process', educationReg, "education")
barPlotJitter('Nudge_type', educationReg, "nudge_education")
table(educationReg$Nudge_type)

#"environment"
environmentReg <- filter(regressionData, Domain=="environment")
barPlotJitter('Cognitive_process', environmentReg, "environment")
environmentReg$Nudge_type  <- factor(environmentReg$Nudge_type , ordered = F)
barPlotJitter('Nudge_type', environmentReg, "nudge_environment")
table(environmentReg$Nudge_type)

#"finance"
financeReg <- filter(regressionData, Domain=="finance")
barPlotJitter('Cognitive_process', financeReg, "finance")
barPlotJitter('Nudge_type', financeReg, "nudge_finance")
table(financeReg$Nudge_type)

#by domain:  "health" 
healthReg <- filter(regressionData, Domain=="health")
barPlotJitter('Cognitive_process', healthReg, "health")
barPlotJitter('Nudge_type', healthReg, "nudge_health")
table(healthReg$Nudge_type)

#"policy"     
policyReg <- filter(regressionData, Domain=="policy")
barPlotJitter('Cognitive_process', policyReg, "policy")
barPlotJitter('Nudge_type', policyReg, "nudge_policy")
table(policyReg$Nudge_type)


####Examine publication bias####
#Create funnel plot
library("metafor")
library("dmetar")

#eggers test
res <- rma(yi = convertedES, sei = se,  data=newMeta)
regtest(res)

#standard funnel plot
res <- rma(yi=convertedES, vi=var, data=newMeta)
funnel(res, ylim=c(0, .25, .5, .75, 1), xlab=expression( paste("Effect size (", italic("d"), ")")))
